Wyjaśnialne uczenie maszynowe – praca domowa 2

Katarzyna Koprowska

In [1]:
import pandas as pd
import numpy as np
import pickle
In [2]:
import warnings
warnings.filterwarnings("ignore")
In [3]:
import matplotlib.pyplot as plt

Wczytanie danych

Wykorzystanym zbirem danych jest Home Equity (HMEQ), zawierający informacje o 5960 klientach banku, którzy otrzymali kredyty hipoteczne.

Na podstawie zbioru próbowałam przewidzieć prawdopodobieństwo defaultu, czyli faktu, że klient będzie zalegał z płatnościami – określa to binarna zmienna BAD (1 oznacza default). Pozostałe 12 zmiennych opisuje m.in. historię kredytową aplikującego, historię zawodową oraz charakterystyki obecnej pożyczki.

Więcej informacji na temat danych można znaleźć pod linkiem https://www.kaggle.com/ajay1735/hmeq-data

In [4]:
hmeq = pd.read_csv("hmeq.csv", error_bad_lines=False)
In [5]:
hmeq_info = {'BAD' : 'client defaulted on loan 0 = loan repaid',
"LOAN" : "Amount of the loan request",
"MORTDUE" : "Amount due on existing mortgage",
"VALUE": "Value of current property",
"REASON": "DebtCon debt consolidation HomeImp = home improvement",
"JOBS" : "occupational categories",
"YOJ": "Years at present job",
"DEROG" : "Number of major derogatory reports",
"DELINQ": "Number of delinquent credit lines",
"CLAGE": "Age of oldest trade line in months",
"NINQ": "Number of recent credit lines",
"CLNO": "Number of credit lines",
"DEBTINC" : "Debt-to-income ratio"}

Przekształcenie danych nienumerycznych na dummy variables

In [6]:
from pandas.api.types import is_numeric_dtype
{column : is_numeric_dtype(hmeq[column]) for column in hmeq.columns}
Out[6]:
{'BAD': True,
 'LOAN': True,
 'MORTDUE': True,
 'VALUE': True,
 'REASON': False,
 'JOB': False,
 'YOJ': True,
 'DEROG': True,
 'DELINQ': True,
 'CLAGE': True,
 'NINQ': True,
 'CLNO': True,
 'DEBTINC': True}
In [7]:
set(hmeq['REASON'])
Out[7]:
{'DebtCon', 'HomeImp', nan}
In [8]:
set(hmeq['JOB'])
Out[8]:
{'Mgr', 'Office', 'Other', 'ProfExe', 'Sales', 'Self', nan}
In [9]:
hmeq = pd.concat([hmeq, pd.get_dummies(hmeq['REASON'], prefix='REASON', dummy_na=True)],axis=1)
hmeq = pd.concat([hmeq, pd.get_dummies(hmeq['JOB'], prefix='JOB', dummy_na=True)],axis=1)
hmeq.drop(['REASON', 'JOB'],axis=1, inplace=True)

Braki danych

In [10]:
hmeq.isna().sum()
Out[10]:
BAD                  0
LOAN                 0
MORTDUE            518
VALUE              112
YOJ                515
DEROG              708
DELINQ             580
CLAGE              308
NINQ               510
CLNO               222
DEBTINC           1267
REASON_DebtCon       0
REASON_HomeImp       0
REASON_nan           0
JOB_Mgr              0
JOB_Office           0
JOB_Other            0
JOB_ProfExe          0
JOB_Sales            0
JOB_Self             0
JOB_nan              0
dtype: int64
In [11]:
hmeq_nonan = hmeq.dropna()
In [12]:
X = hmeq_nonan.iloc[:, 1:]
y = hmeq_nonan.loc[:, "BAD"]
In [13]:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.35, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(X_test, y_test, test_size=0.6, random_state=42)
In [14]:
for data in [X_train, X_test, X_val, y_train,  y_val, y_test]:
    data.reset_index(drop=True, inplace = True)
In [15]:
X_train.shape
Out[15]:
(2284, 20)
In [16]:
metrics = ["accuracy_train", "accuracy_test", "roc_auc_train", "roc_auc_test"]

Model – las losowy

In [17]:
from sklearn.ensemble import RandomForestClassifier
In [18]:
rf_final1 = pickle.load(open("final_nonan_rf.p", "rb"))

Sprawdzenie na zbiorze testowym

In [19]:
from sklearn.metrics import accuracy_score, roc_auc_score
In [20]:
results = {metric : {} for metric in ["accuracy_test", "roc_auc_test"]}
results["accuracy_test"]["RandomForest"] = (accuracy_score(y_test, rf_final1.predict(X_test)))
results["roc_auc_test"]["RandomForest"] = (roc_auc_score(y_test, rf_final1.predict_proba(X_test)[:,1]))
In [21]:
results = pd.DataFrame(results)
In [22]:
results
Out[22]:
accuracy_test roc_auc_test
RandomForest 0.941813 0.895652

Wyjaśnianie

[2. for some selected observation from this dataset, calculate the model predictions for model (1)]

In [23]:
np.random.seed(42)
ind = np.random.randint(len(X_train.index)+1, size=1)[0]
obs = pd.DataFrame(X_train.iloc[ind, :]).T
In [24]:
y_train[obs.index].values
Out[24]:
array([0])
In [25]:
obs.size
Out[25]:
20
In [26]:
rf_final1.predict_proba(obs)
Out[26]:
array([[0.97820614, 0.02179386]])

[3. for an observation selected in (2), calculate the decomposition of model prediction using SHAP, Break Down or both (packages for R: DALEX, iml, packages for python: shap, dalex, piBreakDown).]

Przykładowe wyjaśnienie metodą breakdown.

In [27]:
import dalex
In [28]:
from dalex.explainer import Explainer
In [29]:
exp = Explainer(model = rf_final1, data = X_train, y = rf_final1.predict_proba(X_train)[:,1], model_type = "classification")

Preparation of a new explainer is initiated

  -> label             : not specified, model's class taken instead!
  -> data              : 2284 rows 20 cols
  -> target variable   : 2284 values
  -> predict function  : <function yhat.<locals>.<lambda> at 0x1a18e7d8c8> will be used
  -> predicted values  : min = 0.0, mean = 0.07908669899593004, max = 1.0
  -> residual function : difference between y and yhat
  -> residuals         : min = 0.0, mean = 0.0, max = 0.0
  -> model_info        : package sklearn

A new explainer has been created!
In [30]:
exp_test = Explainer(model = rf_final1, data = X_test, y = rf_final1.predict_proba(X_test)[:,1], model_type = "classification")

Preparation of a new explainer is initiated

  -> label             : not specified, model's class taken instead!
  -> data              : 739 rows 20 cols
  -> target variable   : 739 values
  -> predict function  : <function yhat.<locals>.<lambda> at 0x1a18e7db70> will be used
  -> predicted values  : min = 0.0, mean = 0.07647628756544961, max = 1.0
  -> residual function : difference between y and yhat
  -> residuals         : min = 0.0, mean = 0.0, max = 0.0
  -> model_info        : package sklearn

A new explainer has been created!
In [31]:
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot()

Przykładowe wyjaśnienie metodą SHAP.

In [32]:
import shap
In [33]:
shap.initjs()

X,y = X_train.reset_index(drop=True), y_train.reset_index(drop=True)
model = rf_final1

# explain the model's predictions using SHAP
# (same syntax works for LightGBM, CatBoost, scikit-learn and spark models)
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X, check_additivity=False)

# visualize the first prediction's explanation (use matplotlib=True to avoid Javascript)
shap.force_plot(explainer.expected_value[0], shap_values[0][0,:], X.iloc[0,:])
Out[33]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.
In [34]:
shap.force_plot(explainer.expected_value[0], shap_values[0], X)
Out[34]:
Visualization omitted, Javascript library not loaded!
Have you run `initjs()` in this notebook? If this notebook was from another user you must also trust this notebook (File -> Trust notebook). If you are viewing this notebook on github the Javascript has been stripped for security. If you are using JupyterLab this error is because a JupyterLab extension has not yet been written.

[4. find two observations in the data set, such that they have different most important variables (e.g. age and gender are the most important for observation A, but race and class for observation B)]

In [35]:
important_variables = {}
for i in range(len(X_train.index)):
    obs = pd.DataFrame(X_train.iloc[i, :]).T
    break_down = exp.predict_parts(obs, type = "break_down")
    df = break_down.result.iloc[break_down.result["contribution"].abs().argsort()]
    var_sign = df.loc[:, ["variable_name", "sign"]]
    indexNames = var_sign[var_sign['variable_name'].isin(["intercept", "prediction", ""]) ].index
    var_sign.drop(indexNames , inplace=True)
    important_variables[i] = {"variable_name" : var_sign.loc[:, "variable_name"].values,
                             "sign" : var_sign.loc[:, "sign"].values}
In [36]:
variables = [value['variable_name'] for key, value in important_variables.items()]
In [37]:
signs = [value['sign'] for key, value in important_variables.items()]
In [38]:
variables_2 = [i[-2:] for i in variables]
signs_2 = [i[-2:] for i in signs]

Obserwacjami o różnych najbardziej wpływowych zmiennych okazały się m.in. te o indeksach 2 i 10 ze zbioru treningowego.

In [39]:
ind_1 = 2
ind_2 = 10

Co ciekawe, mają one te same wartości zmiennej zależnej.

In [40]:
print(y_train.reset_index(drop=True)[ind_1])
print(y_train.reset_index(drop=True)[ind_1])

0
0

W przypadku obserwacji o indeksie 2 są to zmienne DEBTINC i DELINQ oznaczające odpowiednio stosunek długu do przychodu i liczbę istniejących kredytów, w których aplikujący zalega ze spłatami, natomiast na wynik obserwacji 10 najbardziej wpłynęły CLAGE i MORTDUE, czyli "wiek" najstarszej linii kredytowej oraz kwota do spłacenia na istniejącym kredycie.

In [41]:
print(variables_2[ind_1])
print(variables_2[ind_2])

['CLAGE' 'CLNO']
['DEBTINC' 'VALUE']
In [42]:
for feature in list(variables_2[ind_1])+list(variables_2[ind_2]):
    print(feature + ": " + hmeq_info[feature])

CLAGE: Age of oldest trade line in months
CLNO: Number of credit lines
DEBTINC: Debt-to-income ratio
VALUE: Value of current property

Możemy przyjrzeć się tym obserwacjom na wykresach. Poniżej pierwsza z nich, czyli obserwacja o indeksie 2.

In [43]:
obs = pd.DataFrame(X_train.iloc[ind_1, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))

Możemy zauważyć, że wartość stosunku długu do przychodu ok. 28 zmniejsza predykcję (czyli zmniejsza prawdopodobieństwo wejścia w default) o 0.055, za to liczba istniejących kredytów równa 1 zwiększa ją o 0.11.

Poniżej wykres rozbicia obserwacji o indeksie 10.

In [44]:
obs = pd.DataFrame(X_train.iloc[ind_2, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))

Widzimy, że najbardziej wpływową zmienną MORTDUE, czyli kwota pozostała do spłacenia na istniejącym kredycie, która będąc równa ok. 50000 zmniejsza predykcję o 0.031. Odwrotny efekt ma zmienna CLAGE (wiek najstarszej linii kredytowej), która przy wartości 93 zwiększa predykcję o 0.026.

[5. select one variable and find two observations in the data set such that for one observation this variable has a positive effect and for the other a negative effect]

In [45]:
vars_signs = [[(k, l) for k, l in zip(i, j)] for i, j in zip(variables, signs)]

Wybraną przez mnie zmienną do porównania jest DEBTINC czyli stosunek kwoty długu do przychodu.

Obserwacjami, które będę porównywać są obserwacje ze zbioru treningowego o indeksach 0 i 1.

In [46]:
ind_1=0
ind_2=1

Mają one inne wartości zmiennej zależnej.

In [47]:
print(y_train[0])
print(y_train[1])

1
0
In [48]:
print(vars_signs[ind_1][-1])
print(vars_signs[ind_2][-1])

('CLAGE', 1.0)
('DELINQ', 1.0)

Wizualizacja rozbicia pierwszej obserwacji.

In [49]:
obs = pd.DataFrame(X_train.iloc[ind_1, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))

Widzimy, że wartość zmiennej DEBTINC równa 34 powoduje obniżenie predykcji o 0.022.

Wizualizacja rozbicia drugiej obserwacji.

In [50]:
obs = pd.DataFrame(X_train.iloc[ind_2, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))

Tutaj z kolei możemy zaobserwować, że wartość tej zmiennej równa 46 wpływa na znaczne podwyższenie predykcji.

[6. train a second model (of any class, neural nets, linear, other boosting) and find an observation for which BD/shap attributions are different between the models]

In [51]:
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(solver="liblinear", random_state=42, intercept_scaling= 1.3).fit(X_train, y_train)
In [52]:
exp_lr = Explainer(model = lr, data = X_train, y = lr.predict_proba(X_train)[:,1], model_type = "classification")

Preparation of a new explainer is initiated

  -> label             : not specified, model's class taken instead!
  -> data              : 2284 rows 20 cols
  -> target variable   : 2284 values
  -> predict function  : <function yhat.<locals>.<lambda> at 0x1a311159d8> will be used
  -> predicted values  : min = 3.367270657199028e-06, mean = 0.09338772022355352, max = 0.44416952305359814
  -> residual function : difference between y and yhat
  -> residuals         : min = 0.0, mean = 0.0, max = 0.0
  -> model_info        : package sklearn

A new explainer has been created!
In [ ]:
important_variables_lr = {}
for i in range(len(X_train.index)):
    obs = pd.DataFrame(X_train.iloc[i, :]).T
    break_down = exp_lr.predict_parts(obs, type = "break_down")
    df = break_down.result.iloc[break_down.result["contribution"].abs().argsort()]
    var_sign = df.loc[:, ["variable_name", "sign"]]
    indexNames = var_sign[var_sign['variable_name'].isin(["intercept", "prediction", ""]) ].index
    var_sign.drop(indexNames , inplace=True)
    important_variables_lr[i] = {"variable_name" : var_sign.loc[:, "variable_name"].values,
                             "sign" : var_sign.loc[:, "sign"].values}
In [ ]:
variables_lr = [value['variable_name'] for key, value in important_variables_lr.items()]
In [ ]:
signs_lr = [value['sign'] for key, value in important_variables_lr.items()]
In [ ]:
variables_2_lr = [i[-2:] for i in variables_lr]
signs_2_lr = [i[-2:] for i in signs_lr]
In [ ]:
ind = 0
In [ ]:
obs = pd.DataFrame(X_train.iloc[ind, :]).T
break_down = exp.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))
In [ ]:
obs = pd.DataFrame(X_train.iloc[ind, :]).T
break_down = exp_lr.predict_parts(obs, type = "break_down")
break_down.plot(max_vars=len(X_train.columns))

W przypadku lasów losowych najbardziej wpływową zmienną jest znany nam już DEBTINC, czyli stosunek długu do dochodu, który znacząco obniża predykcję.

W regresji logistycznej ma on dużo mniejszy wpływ (co więcej – w drugą stronę), za to na prowadzenie wychodzi zmienna YOJ oznaczająca liczbę lat w obecnej pracy, która w modelu lasów losowych miała bardzo niewielki wpływ, którego zwrot także różnił się od zwrotu w modelu regresji.

Podsumowanie

Jak zostało pokazane, wyjaśnienia poszczególnych obserwacji w obrębie modelu mogą być bardzo różne, nawet gdy wartości zmiennej zależnej są takie same. Można też było zauważyć, że choć dana wartość zmiennej obniża predykcję, niewielkie, zdawałoby się, zwiększenie tej wartości skutkuje całkowicie różnym wpływem na predykcję.

Dodatkowo, różnice mogą występować także w obrębie jednej obserwacji -- w zależności od modelu zmienne wpływające na jej predykcję mogą się bardzo różnić.